      MODULE GWFSWIMODULE
        TYPE TSWIDE4
          INTEGER  ::MXITER,NODES,NHALFU,NHALFL,NBWGRD
          INTEGER  ::MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR
          INTEGER  ::NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM
          INTEGER  ::NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
          REAL     ::ACCLDE4,HCLOSEDE4,DELTL
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   ::IUPPNT
          INTEGER,          ALLOCATABLE, DIMENSION(:,:,:) ::IEQPNT
          REAL,             ALLOCATABLE, DIMENSION(:,:)   ::AU
          REAL,             ALLOCATABLE, DIMENSION(:,:)   ::AL
          REAL,             ALLOCATABLE, DIMENSION(:)     ::D4B
          REAL,             ALLOCATABLE, DIMENSION(:)     ::HDCGDE4
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   ::LRCHDE4
        END TYPE TSWIDE4
        TYPE TSWIPCG
          INTEGER ::MXITER,NODES
          INTEGER ::ITER1,NPCOND,NBPOL,IPRPCG,MUTPCG,NITER
          REAL    ::HCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG
          REAL    ::DAMPPCGT
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) ::VPCG
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) ::SS
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) ::P
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) ::HPCG
          REAL,             ALLOCATABLE, DIMENSION(:,:,:) ::CD
          REAL,             ALLOCATABLE, DIMENSION(:,:,:) ::HCSV
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   ::LHCH
          REAL,             ALLOCATABLE, DIMENSION(:)     ::HCHG
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   ::LRCHPCG
          REAL,             ALLOCATABLE, DIMENSION(:)     ::RCHG
          INTEGER,          ALLOCATABLE, DIMENSION(:)     ::IT1
        END TYPE TSWIPCG
        INTEGER, SAVE, POINTER :: NPLN,ISTRAT,ISWIZT,NPRS,NZONES
        INTEGER, SAVE, POINTER :: NLAYSWI
!        INTEGER, SAVE, POINTER :: NODES,NHALFU,NHALFL
        REAL, SAVE, POINTER    :: USEDE45
!        INTEGER, SAVE, POINTER :: ID4DIM,ID4DIR,NBWGRD
!        INTEGER, SAVE, POINTER :: ITMX,MXITER,MXUP,MXLOW,MXBW,MXEQ
!        INTEGER, SAVE, POINTER :: MUTD4,NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
!        REAL, SAVE, POINTER    :: DELTL
        REAL, SAVE, POINTER    :: TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA
!        REAL, SAVE, DIMENSION(:,:), POINTER :: AU
!        INTEGER, SAVE, DIMENSION(:,:), POINTER :: IUPP
!        REAL, SAVE, DIMENSION(:,:), POINTER :: AL
!        INTEGER, SAVE, DIMENSION(:,:,:), POINTER :: IEQP
!        REAL, SAVE, DIMENSION(:), POINTER :: SWID4B
!        INTEGER, SAVE, DIMENSION(:,:), POINTER :: SWILRCH
!        REAL, SAVE, DIMENSION(:), POINTER :: SWIHDCG
        INTEGER, SAVE, DIMENSION(:,:), POINTER :: IBO
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWIHCOF
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWISOLVCR
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWISOLVCC
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWISOLVCV
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: ZETA
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: ZETAOLD
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: SSZ
        REAL, SAVE, DIMENSION(:), POINTER :: EPS
        REAL, SAVE, DIMENSION(:), POINTER :: NUS
        REAL, SAVE, DIMENSION(:), POINTER :: DELNUS
        REAL, SAVE, DIMENSION(:), POINTER :: NUPLANE
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICR
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICC
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICUMCR
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICUMCC
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: NUTOP
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: NUBOT
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: BRHS
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: QZEXTRA
        DOUBLE PRECISION, SAVE, DIMENSION(:,:), POINTER :: DUM
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: RHSFRESH
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: HCOFOLD
        INTEGER, SAVE, DIMENSION(:,:,:,:), POINTER :: IPLPOS
        INTEGER, SAVE, DIMENSION(:,:,:), POINTER :: IZONENR
        REAL,    SAVE,   POINTER, DIMENSION(:,:,:) ::SC1
C---------STORAGE FOR SOLVERS
        TYPE (TSWIDE4), SAVE, POINTER :: SWIDE4        
        TYPE (TSWIPCG), SAVE, POINTER :: SWIPCG        

        TYPE GWFSWITYPE
          INTEGER, POINTER :: NPLN,ISTRAT,ISWIZT,NPRS,NZONES
          INTEGER, POINTER :: NLAYSWI
!          INTEGER, POINTER :: NODES,NHALFU,NHALFL
          REAL, POINTER    :: USEDE45
!          INTEGER, POINTER :: ID4DIM,ID4DIR,NBWGRD
!          INTEGER, POINTER :: ITMX,MXITER,MXUP,MXLOW,MXBW,MXEQ
!          INTEGER, POINTER :: MUTD4,NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
!          REAL, POINTER    :: DELTL
          REAL, POINTER    :: TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA
!          REAL, DIMENSION(:,:), POINTER :: AU
!          INTEGER, DIMENSION(:,:), POINTER :: IUPP
!          REAL, DIMENSION(:,:), POINTER :: AL
!          INTEGER, DIMENSION(:,:,:), POINTER :: IEQP
!          REAL, DIMENSION(:), POINTER :: SWID4B
!          INTEGER, DIMENSION(:,:), POINTER :: SWILRCH
!          REAL, DIMENSION(:), POINTER :: SWIHDCG
          INTEGER, DIMENSION(:,:), POINTER :: IBO
          REAL, DIMENSION(:,:), POINTER :: SWIHCOF
          REAL, DIMENSION(:,:), POINTER :: SWISOLVCR
          REAL, DIMENSION(:,:), POINTER :: SWISOLVCC
          REAL, DIMENSION(:,:), POINTER :: SWISOLVCV
          REAL, DIMENSION(:,:,:,:), POINTER :: ZETA
          REAL, DIMENSION(:,:,:,:), POINTER :: ZETAOLD
          REAL, DIMENSION(:,:,:), POINTER :: SSZ
          REAL, DIMENSION(:), POINTER :: EPS
          REAL, DIMENSION(:), POINTER :: NUS
          REAL, DIMENSION(:), POINTER :: DELNUS
          REAL, DIMENSION(:), POINTER :: NUPLANE
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICR
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICC
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICUMCR
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICUMCC
          REAL, DIMENSION(:,:,:), POINTER :: NUTOP
          REAL, DIMENSION(:,:,:), POINTER :: NUBOT
          REAL, DIMENSION(:,:,:), POINTER :: BRHS
          REAL, DIMENSION(:,:,:), POINTER :: QZEXTRA
          DOUBLE PRECISION, DIMENSION(:,:), POINTER :: DUM
          REAL, DIMENSION(:,:,:), POINTER :: RHSFRESH
          REAL, DIMENSION(:,:,:), POINTER :: HCOFOLD
          INTEGER, DIMENSION(:,:,:,:), POINTER :: IPLPOS
          INTEGER, DIMENSION(:,:,:), POINTER :: IZONENR
          REAL, DIMENSION(:,:,:), POINTER ::SC1
C-----------STORAGE FOR SOLVERS
          TYPE (TSWIDE4), POINTER :: SWIDE4        
          TYPE (TSWIPCG), POINTER :: SWIPCG        
        END TYPE
        TYPE(GWFSWITYPE), SAVE:: GWFSWIDAT(10)
      END MODULE GWFSWIMODULE


      SUBROUTINE GWF2SWI1AR(In,Ibcf,Ilpf,Ihuf,Igrid)
C-----VERSION 3 26SEP2010 GWF2SWI1AR
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR SEA WATER INTRUSION PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                      LBOTM,BOTM  
      USE GWFBCFMODULE,ONLY:SC1B=>SC1
      USE GWFLPFMODULE,ONLY:SC1L=>SC1
      USE GWFHUFMODULE,ONLY:SC1H=>SC1
      USE GWFSWIMODULE
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: In
      INTEGER, INTENT(IN) :: Ibcf
      INTEGER, INTENT(IN) :: Ilpf
      INTEGER, INTENT(IN) :: Ihuf
      INTEGER, INTENT(IN) :: Igrid
C     + + + LOCAL DEFINITIONS + + +
      INTEGER :: i, j, k
      INTEGER :: iz, kk
      INTEGER :: itmem
      REAL :: d
      REAL :: bbot, ttop, z
	CHARACTER*24 ZETANAME
C     ------------------------------------------------------------------
C     ------------------------------------------------------------------
C
C
        ALLOCATE(NPLN,ISTRAT,ISWIZT,NPRS,NZONES)
        ALLOCATE(NLAYSWI)
!        ALLOCATE(NODES,NHALFU,NHALFL)
        ALLOCATE(USEDE45)
!        ALLOCATE(ID4DIM,ID4DIR,NBWGRD)
!        ALLOCATE(ITMX,MXITER,MXUP,MXLOW,MXBW,MXEQ)
!        ALLOCATE(MUTD4,NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW)
!        ALLOCATE(DELTL)

        ALLOCATE(TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA)


C-------IDENTIFY PACKAGE AND INITIALIZE
      WRITE(IOUT,1) In
    1 FORMAT(1X,/1X,'SWI1 -- SWI PACKAGE, VERSION 3.1, 25-Aug-2004',
     1' INPUT READ FROM UNIT',I3)
C
C-------READ & PRINT NUMBER OF PLANES
C------- & DENSITY OPTIONS AND UNIT OR FLAG FOR SAVING ZETA
      READ(In,2) NPLN,ISTRAT,ISWIZT,NPRS
    2 FORMAT(4I10)
C
C-------FLAG FOR USE OF PCG SOLVER INSTEAD OF DE4
      USEDE45 = 1
      IF ( NPLN.LT.0 ) THEN
        USEDE45 = 0
        NPLN = ABS( NPLN )
      END IF
      IF ( USEDE45.EQ.1 ) THEN
        WRITE (IOUT,'(/1XA//)') 
     &    'DE45 SOLVER USED TO SOLVE FOR SWI ZETA SURFACES'
      ELSE IF ( USEDE45.EQ.0 ) THEN
        WRITE (IOUT,'(/1XA//)') 
     &    'PCG2 SOLVER USED TO SOLVE FOR SWI ZETA SURFACES'
      END IF

      IF (ISWIZT.GT.0) WRITE(IOUT,3) ISWIZT
    3 FORMAT(1X,'ZETA WILL BE SAVED ON UNIT',I3)
C
C3-----CHECK TO SEE THAT NUMBER OF PLANES IS AT LEAST 1,
C3-----PRINT VALUE
      IF(NPLN.GT.0) THEN
        WRITE(IOUT,6) NPLN
      ELSE
        WRITE (IOUT,7)
        CALL USTOP('ABORTING, NUMBER OF PLANES LESS THAN 1...')
      ENDIF 
    6 FORMAT(1X,'TOTAL NUMBER OF PLANES: ',I3)
    7 FORMAT(1X,'ABORTING, NUMBER OF PLANES LESS THAN 1...')
C
C4-----WRITE DENSITY OPTIONS TO LIST FILE
	IF(ISTRAT.EQ.1) WRITE(IOUT,8)
    8 FORMAT(1X,'ISTRAT 1 -- STRATIFIED FLOW')
	IF(ISTRAT.EQ.0) WRITE(IOUT,9)
    9 FORMAT(1X,'ISTRAT 0 -- VARIABLE DENSITY FLOW')

      NZONES=NPLN+1;
      NLAYSWI=1
C
C5------COPIED FROM DE45 PACKAGE SUBSTITUTE NLAY BY 1
C5------CALCULATE BAND WIDTH
!      NODES=NCOL*NROW*NLAYSWI
!      NHALFU=(NODES-1)/2 + 1
!      NHALFL=NODES-NHALFU
!      
!      ID4DIM=7
!      !MODFLOW-2005 - START
!      DELTL = 0.
!      NBWL = 0
!      NUPL = 0
!      NLOWL = 0
!      !MODFLOW-2005 - END
!      IF(NLAYSWI.LE.NCOL .AND. NLAYSWI.LE.NROW) THEN
!         IF(NLAYSWI.EQ.1) ID4DIM=5
!         IF(NCOL.GE.NROW) THEN
!            ID4DIR=1
!            NBWGRD=NROW*NLAYSWI+1
!         ELSE
!            ID4DIR=2
!            NBWGRD=NCOL*NLAYSWI+1
!         END IF
!      ELSE IF(NROW.LE.NCOL .AND. NROW.LE.NLAYSWI) THEN
!         IF(NROW.EQ.1) ID4DIM=5
!         IF(NCOL.GE.NLAYSWI) THEN
!            ID4DIR=3
!            NBWGRD=NROW*NLAYSWI+1
!         ELSE
!            ID4DIR=4
!            NBWGRD=NROW*NCOL+1
!         END IF
!      ELSE
!         IF(NCOL.EQ.1) ID4DIM=5
!         IF(NROW.GE.NLAYSWI) THEN
!            ID4DIR=5
!            NBWGRD=NCOL*NLAYSWI+1
!         ELSE
!            ID4DIR=6
!            NBWGRD=NCOL*NROW+1
!         END IF
!      END IF
!      MXUP=NHALFU
!      MXLOW=NHALFL
!      !SWI VERSION 2
!	!MXBW=NBWGRD + 4
!      !SWI VERSION 2
!
!      !MODFLOW-2005 SWI VERSION 3
!	MXBW=NBWGRD
!      !MODFLOW-2005 SWI VERSION 3
!      MXEQ=MXUP+MXLOW
C
      ALLOCATE(ZETA(NCOL,NROW,NLAY,NZONES+1))
      ALLOCATE(SSZ(NCOL,NROW,NLAY))
      ALLOCATE(IZONENR(NCOL,NROW,NLAY))
      ALLOCATE(NUS(NZONES))
	ALLOCATE(NUPLANE(NZONES+1))
	ALLOCATE(DELNUS(NZONES))
	ALLOCATE(EPS(NZONES))

	ALLOCATE(IPLPOS(NCOL,NROW,NLAY,NZONES))
      ALLOCATE(RHSFRESH(NCOL,NROW,NLAY))
	ALLOCATE(SWICR(NCOL,NROW,NLAY,NZONES))
	ALLOCATE(SWICC(NCOL,NROW,NLAY,NZONES))
      ALLOCATE(SWICUMCR(NCOL,NROW,NLAY,NZONES))
	ALLOCATE(SWICUMCC(NCOL,NROW,NLAY,NZONES))
	ALLOCATE(NUTOP(NCOL,NROW,NLAY),NUBOT(NCOL,NROW,NLAY))
	ALLOCATE(QZEXTRA(NCOL,NROW,NLAY))

      ALLOCATE(BRHS(NCOL,NROW,NZONES))
	ALLOCATE(HCOFOLD(NCOL,NROW,NLAY))
	ALLOCATE(ZETAOLD(NCOL,NROW,NLAY,NZONES+1))
      ALLOCATE(IBO(NCOL,NROW),SWIHCOF(NCOL,NROW))
	ALLOCATE(SWISOLVCR(NCOL,NROW))
	ALLOCATE(SWISOLVCC(NCOL,NROW))
	ALLOCATE(SWISOLVCV(NCOL,NROW))
	ALLOCATE(DUM(NCOL,NROW))
	
!	ALLOCATE(AU(ID4DIM,MXUP))
!	ALLOCATE(AL(MXBW,MXLOW))
!	ALLOCATE(IEQP(NCOL,NROW,NLAYSWI))
!	ALLOCATE(IUPP(ID4DIM,MXUP))
!	ALLOCATE(SWID4B(MXEQ))
!	ALLOCATE(SWILRCH(3,ITMX))
!	ALLOCATE(SWIHDCG(ITMX))
C
C-------ALLOCATE STORAGE FOR SOLVER DATA
      ALLOCATE(SWIDE4,SWIPCG)
C-------ALLOCATE SPACE FOR THE DE4 ARRAYS
      IF ( USEDE45.EQ.1 ) THEN
        SWIDE4%ITMX = 1
!---------SPECIFY DEFAULT PARAMETERS - MODIFY LATER TO ALLOW USERS TO SPECIFY IN SWI FILE 
        SWIDE4%MXITER = SWIDE4%ITMX
	  SWIDE4%NITERDE4 = 1
	  SWIDE4%ACCLDE4 = 1.0
	  SWIDE4%HCLOSEDE4 = 1.0E-05
	  SWIDE4%IPRD4 = 0
	  SWIDE4%IFREQ = 3
	  SWIDE4%DELTL = 1.
C         INITIALIZE DIMENSION PARAMETERS
        SWIDE4%NODES=NCOL*NROW*NLAYSWI
        SWIDE4%NHALFU=(SWIDE4%NODES-1)/2 + 1
        SWIDE4%NHALFL=SWIDE4%NODES-SWIDE4%NHALFU
        SWIDE4%NBWL = 0
        SWIDE4%NUPL = 0
        SWIDE4%NLOWL = 0
        SWIDE4%ID4DIM = 7
C         CALCULATE SOLVER DIMENSIONS
        IF(NLAYSWI.LE.NCOL .AND. NLAYSWI.LE.NROW) THEN
           IF(NLAYSWI.EQ.1) SWIDE4%ID4DIM=5
           IF(NCOL.GE.NROW) THEN
              SWIDE4%ID4DIR=1
              SWIDE4%NBWGRD=NROW*NLAYSWI+1
           ELSE
              SWIDE4%ID4DIR=2
              SWIDE4%NBWGRD=NCOL*NLAYSWI+1
           END IF
        ELSE IF(NROW.LE.NCOL .AND. NROW.LE.NLAYSWI) THEN
           IF(NROW.EQ.1) SWIDE4%ID4DIM=5
           IF(NCOL.GE.NLAYSWI) THEN
              SWIDE4%ID4DIR=3
              SWIDE4%NBWGRD=NROW*NLAYSWI+1
           ELSE
              SWIDE4%ID4DIR=4
              SWIDE4%NBWGRD=NROW*NCOL+1
           END IF
        ELSE
           IF(NCOL.EQ.1) SWIDE4%ID4DIM=5
           IF(NROW.GE.NLAYSWI) THEN
              SWIDE4%ID4DIR=5
              SWIDE4%NBWGRD=NCOL*NLAYSWI+1
           ELSE
              SWIDE4%ID4DIR=6
              SWIDE4%NBWGRD=NCOL*NROW+1
           END IF
        END IF
        SWIDE4%MXUP=SWIDE4%NHALFU
        SWIDE4%MXLOW=SWIDE4%NHALFL
	  SWIDE4%MXBW=SWIDE4%NBWGRD
        SWIDE4%MXEQ=SWIDE4%MXUP+SWIDE4%MXLOW
C         ALLOCATE ARRAYS      
        ALLOCATE (SWIDE4%AU(SWIDE4%ID4DIM,SWIDE4%MXUP))
        ALLOCATE (SWIDE4%IUPPNT(SWIDE4%ID4DIM,SWIDE4%MXUP))
        ALLOCATE (SWIDE4%AL(SWIDE4%MXBW,SWIDE4%MXLOW))
        ALLOCATE (SWIDE4%IEQPNT(NCOL,NROW,NLAYSWI))
        ALLOCATE (SWIDE4%D4B(SWIDE4%MXEQ))
        ALLOCATE (SWIDE4%LRCHDE4(3,SWIDE4%ITMX))
        ALLOCATE (SWIDE4%HDCGDE4(SWIDE4%ITMX))
C        INITIALIZE FLOAT ARRAYS
        SWIDE4%AU = 0.0
        SWIDE4%AL = 0.0
        SWIDE4%D4B = 0.0
      ELSE
!---------SPECIFY DEFAULT PARAMETERS - MODIFY LATER TO ALLOW USERS TO SPECIFY IN SWI FILE 
        SWIPCG%MXITER = 50
        SWIPCG%ITER1 = 25
        SWIPCG%NPCOND = 1  
        SWIPCG%HCLOSEPCG = 1.0E-05
        SWIPCG%RCLOSEPCG = 1.0E-05
        SWIPCG%RELAXPCG = 1.0
        SWIPCG%NBPOL = 0
        SWIPCG%IPRPCG = 999
        SWIPCG%MUTPCG = 3
        SWIPCG%DAMPPCG = 1.0
        SWIPCG%DAMPPCGT = 1.0
C         INITIALIZE DIMENSION PARAMETERS
        SWIPCG%NODES=NCOL*NROW*NLAYSWI
C-------ALLOCATE SPACE FOR THE PCG ARRAYS
        ALLOCATE (SWIPCG%VPCG(NCOL,NROW,NLAYSWI))
        ALLOCATE (SWIPCG%SS(NCOL,NROW,NLAYSWI))
        ALLOCATE (SWIPCG%P(NCOL,NROW,NLAYSWI))
        ALLOCATE (SWIPCG%HPCG(NCOL,NROW,NLAYSWI))
        ALLOCATE (SWIPCG%CD(NCOL,NROW,NLAYSWI))
        IF(SWIPCG%NPCOND.EQ.2) THEN
           ALLOCATE (SWIPCG%HCSV(NCOL,NROW,NLAYSWI))
        ELSE
           ALLOCATE (SWIPCG%HCSV(1,1,1))
        END IF
        itmem=SWIPCG%MXITER*SWIPCG%ITER1
        ALLOCATE (SWIPCG%HCHG(itmem))
        ALLOCATE (SWIPCG%LHCH(3,itmem))
        ALLOCATE (SWIPCG%RCHG(itmem))
        ALLOCATE (SWIPCG%LRCHPCG(3,itmem))
        ALLOCATE (SWIPCG%IT1(itmem))
C         INITIALIZE ARRAYS
        SWIPCG%VPCG = 0.0        
        SWIPCG%SS = 0.0
        SWIPCG%P = 0.0
        SWIPCG%HPCG = 0.0
        SWIPCG%CD = 0.0
        SWIPCG%HCSV = 0.0
      END IF
C
C-------POINT SWI SC1 TO APPROPRIATE STORAGE PARAMETER
      IF ( Ibcf.GT.0 ) THEN
        SC1 => SC1B
      ELSE IF ( Ilpf.GT.0 ) THEN
        SC1 => SC1L
      ELSE IF ( Ihuf.GT.0 ) THEN
        SC1 => SC1H
      ELSE
        WRITE (IOUT,2100)
        CALL USTOP('FLOW PACKAGE SPECIFIED INCONSISTENT WITH SWI') 
      END IF
2100  FORMAT(//1X,'SWI PROCESS REQUIRES USE OF THE BCF, LPF,'
     2       1X,'OR HUF FLOW PACKAGES',//)
C
C-------READ SWI DATA - MOVED FROM SWI1RP
C
C-------READ SIMULATION PARAMETERS
      READ(IN,'(F10.0,F10.0,F10.0,F10.0)') 
	1     TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA


	IF (ISTRAT.EQ.1) THEN
C
C-------READ NU FOR EACH ZONE (NZONES)
       CALL U1DREL(NUS,'                  NUZONE',NZONES,IN,IOUT)

C-------SET EPS EQUAL TO ZERO
	 DO iz=1,NZONES
	  EPS(iz)=0.0
	 END DO

	ELSEIF (ISTRAT.EQ.0) THEN
C
C-------READ NUPLANE FOR EACH PLANE (NZONES+1)
       CALL U1DREL(NUPLANE,'                 NUPLANE',NZONES+1,IN,IOUT)
C
C-------CALCULATE NU AND EPS FROM NUPLANE
	 DO iz=1,NZONES
	  NUS(iz)=0.5*(NUPLANE(iz)+NUPLANE(iz+1))
	  EPS(iz)=(NUPLANE(iz+1)-NUPLANE(iz))/6
	 END DO

	END IF
C	
C-------CALCULATE DELNUS FROM NUS
	DELNUS(1)=NUS(1);
	DO iz=2,NZONES
	 DELNUS(iz)=(NUS(iz)-NUS(iz-1))
	END DO
C
C-------READ ZETA FOR EACH PLANE	
	DO 55 iz=2,NZONES
       DO 55 k=1,NLAY
		  kk=k
	     WRITE(ZETANAME,111) iz-1
  111      FORMAT('           ZETA PLANE ',I2)

            CALL U2DREL(ZETA(1,1,k,iz),ZETANAME,
	1	  NROW,NCOL,kk,IN,IOUT)

C
C-------CHECK SOMETHING
	d=0.001
	DO  j=1,NCOL
       DO  i=1,NROW
		  bbot=BOTM(j,i,LBOTM(k))
          ttop=BOTM(j,i,LBOTM(k)-1)	            
			z=ZETA(j,i,k,iz)
			IF(z.NE.bbot) THEN		
				IF(z.LT.(bbot+d)) THEN
					ZETA(j,i,k,iz)=bbot
				ENDIF
			ENDIF
			IF(z.GT.(ttop-d)) ZETA(j,i,k,iz)=ttop

	 ENDDO
	ENDDO

   55 CONTINUE
C	 
C-------READ SSZ FOR EACH LAYER	
      DO 112 k=1,NLAY
		  kk=k
            CALL U2DREL(SSZ(1,1,k),'                     SSZ',
	1	  NROW,NCOL,kk,IN,IOUT)
  112 CONTINUE
C
C-------READ IZONENR FOR EACH LAYER	

      DO 113 k=1,NLAY
        kk=k
      CALL U2DINT(IZONENR(1,1,kk), '                 IZONENR',
     1	NROW,NCOL,kk,IN,IOUT)
  113 CONTINUE
C
C-------SET POINTERS FOR GRID
      CALL SGWF2SWI1PSV(Igrid)
C
C-------RETURN
      RETURN
      END SUBROUTINE GWF2SWI1AR

C
C-------
	SUBROUTINE GWFSWI1HC(Igrid)
C-----VERSION 3 26SEP2010 GWFSWI1HC
C     ******************************************************************
C     COPY HCOF TO SWI HCOFOLD PRIOR TO ADDING MODFLOW BOUNDARY 
C     CONDITION TERMS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:HCOF
      USE GWFSWIMODULE,ONLY:HCOFOLD
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: Igrid
C     + + + LOCAL DEFINITIONS + + +
C     ------------------------------------------------------------------
C
C-------SET SWI POINTERS FOR GRID
      CALL SGWF2SWI1PNT(Igrid)
C
C-------USE WHOLE ARRAY OPERATION TO SET HCOFOLD TO HCOF
	HCOFOLD = HCOF
C-------RETURN
      RETURN
      END SUBROUTINE GWFSWI1HC

!      SUBROUTINE SWI1FM(BOTM,NBOTM,EPS,NUS,DELNUS,IBOUND,HNEW,CR,CC,
!     &  CV,HCOF,RHS,NCOL,NROW,NLAY,NZONES,IPLPOS,ZETA,SWICR,SWICC,
!     &  SWICUMCR,SWICUMCC,DELR,DELC,NUTOP,NUBOT,QZEXTRA,RHSFRESH,IOUT)

      SUBROUTINE GWF2SWI1FM(Igrid)
C-----VERSION 3 26SEP2010 GWF2SWI1FM
C     ******************************************************************
C     ADD SWI TERMS TO RHS AND HCOF
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                      LBOTM,BOTM,
     3                      CR,CC,CV,HCOF,RHS,
     4                      DELR,DELC,IBOUND,HNEW  
      USE GWFSWIMODULE
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: Igrid
C     + + + LOCAL DEFINITIONS + + +
      INTEGER :: i, j, k
      INTEGER :: iz
      INTEGER :: izrev
      REAL :: bbot, ttop, hhd
      REAL :: sumthickrf, sumthickff
      REAL :: dzeta1, dzeta2
      REAL :: thickrf, thickff
      REAL :: sumswicr, sumswicc
      REAL :: headdiff

C     ------------------------------------------------------------------
C
C-------SET POINTERS FOR GRID
      CALL SGWF2SWI1PNT(Igrid)
C
C--------SET FIRST AND LAST PLANE EQUAL TO TOP AND BOTTOM OF THE AQUIFER
      DO 90 k=1,NLAY
       DO 90 i=1,NROW
        DO 90 j=1,NCOL
         bbot=BOTM(j,i,LBOTM(k))
         ttop=BOTM(j,i,LBOTM(k)-1)
         IF(bbot.GT.ttop) THEN
            WRITE(IOUT,35) k,i,j
   35       FORMAT(1X,/1X,'Negative cell thickness at (layer,row,col)',
     1      I4,',',I4,',',I4)
            WRITE(IOUT,36) ttop,bbot
   36       FORMAT(1X,'Top elevation, bottom elevation:',1P,2G13.5)
            STOP
         END IF
!         IF ((LAYCON(k).NE.0).OR.(LAYCON(k).NE.0)) THEN
!            HHD=HNEW(j,i,k)
!            IF(HHD.LT.ttop) ttop=HHD
!         END IF
         hhd=HNEW(j,i,k)
         IF(hhd.LT.ttop) ttop=hhd
         ZETA(j,i,k,1)=ttop
	   ZETA(j,i,k,NZONES+1)=bbot
   90  CONTINUE
C
C-----COMPUTE HORIZONTAL (CUMULATIVE) CONDUCTANCE FOR EACH ZONE
      DO 120 k=1,NLAY
       DO 120 i=1,NROW
        DO 120 j=1,NCOL
   	    sumthickrf=0
          sumthickff=0
         DO iz=1,NZONES
 	    IF(j.NE.NCOL) THEN
	     dzeta1=ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1);
	     dzeta2=ZETA(j+1,i,k,iz)-ZETA(j+1,i,k,iz+1);
		 IF (dzeta1.LE.0 .OR. dzeta2.LE.0) THEN
            thickrf=0
	     ELSE
	      thickrf=(dzeta1*DELR(j+1)+dzeta2*DELR(j))/
     &	                  (DELR(j)+DELR(j+1))
	     END IF
	     SWICR(j,i,k,iz)=CR(j,i,k)*thickrf
		 sumthickrf=sumthickrf+thickrf
	    ELSE
	     SWICR(j,i,k,iz)=0
           sumthickrf=1
	    END IF
		IF(i.NE.NROW) THEN
		 dzeta1=ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)
           dzeta2=ZETA(j,i+1,k,iz)-ZETA(j,i+1,k,iz+1)
		 IF (dzeta1.LE.0 .OR. dzeta2.LE.0) THEN
            thickff=0
	     ELSE 
	      thickff=(dzeta1*DELC(i+1)+dzeta2*DELC(i))/
     &	                  (DELC(i)+DELC(i+1))
	     END IF
	     SWICC(j,i,k,iz)=CC(j,i,k)*thickff
	     sumthickff=sumthickff+thickff
	    ELSE
           SWICC(j,i,k,iz)=0
	     sumthickff=1
	    END IF
	   END DO 
         DO iz=1,NZONES
	   
CFLEVO	
		IF (sumthickrf.NE.0) THEN
C			sumthickrf=1
	    SWICR(j,i,k,iz)=SWICR(j,i,k,iz)/sumthickrf
		ENDIF
		
		IF (sumthickff.NE.0) THEN
C			sumthickrf=1
			SWICC(j,i,k,iz)=SWICC(j,i,k,iz)/sumthickff
		ENDIF
		
CFLEVO				
	    END DO
	    sumswicr=0
          sumswicc=0
         DO iz=1,NZONES    
          sumswicr=sumswicr+SWICR(j,i,k,NZONES+1-iz)
		  sumswicc=sumswicc+SWICC(j,i,k,NZONES+1-iz)
          SWICUMCR(j,i,k,NZONES+1-iz)=sumswicr
	    SWICUMCC(j,i,k,NZONES+1-iz)=sumswicc
         END DO
  120 CONTINUE
C
C-------SET IPLPOS (1=top, 2=bottom, 0=in between)
      DO 130 k=1,NLAY
       DO 130 i=1,NROW
        DO 130 j=1,NCOL
	   IPLPOS(j,i,k,1)=0
         DO 130 iz=2,NZONES
	    IF (ZETA(j,i,k,iz).GE.ZETA(j,i,k,1)) THEN
	     IPLPOS(j,i,k,iz)=1
          ELSEIF (ZETA(j,i,k,iz).LE.ZETA(j,i,k,NZONES+1)) THEN
           IPLPOS(j,i,k,iz)=2
	    ELSE
	     IPLPOS(j,i,k,iz)=0
	    END IF
  130 CONTINUE
C
C-------SET ZONE CONDUCTANCES TO ZERO IF NOT CELL ON EITHER SIDE HAS ACTIVE ZONE
      DO 140 k=1,NLAY
       DO 140 i=1,NROW
        DO 140 j=1,NCOL
         DO 140 iz=2,NZONES-1
	    IF (j.NE.NCOL) THEN
           IF ((IPLPOS(j,i,k,iz)+IPLPOS(j+1,i,k,iz)+
	1	        IPLPOS(j,i,k,iz+1)+IPLPOS(j+1,i,k,iz+1)).NE.0) THEN
            SWICR(j,i,k,iz)=0
	     END IF
		END IF 
	    IF (i.NE.NROW) THEN
           IF ((IPLPOS(j,i,k,iz)+IPLPOS(j,i+1,k,iz)+
	1	        IPLPOS(j,i,k,iz+1)+IPLPOS(j,i+1,k,iz+1)).NE.0) THEN
            SWICC(j,i,k,iz)=0
	     END IF
	    END IF
  140 CONTINUE
C
C--------COMPUTE NUTOP AND NUBOT
      DO 150 k=1,NLAY
       DO 150 i=1,NROW
        DO 150 j=1,NCOL
	   NUTOP(j,i,k)=NUS(1)
         NUBOT(j,i,k)=NUS(NZONES)
	   DO 150 iz=2,NZONES
	    IF (IPLPOS(j,i,k,iz).EQ.1) THEN
	     NUTOP(j,i,k)=NUTOP(j,i,k)+DELNUS(iz)
	    END IF
	     izrev=NZONES-iz+2
	    IF (IPLPOS(j,i,k,izrev).EQ.2) THEN
	     NUBOT(j,i,k)=NUBOT(j,i,k)-DELNUS(izrev)
	    END IF
  150 CONTINUE
C     
C--------COMPUTE QZEXTRA
	DO 160 k=1,NLAY
       DO 160 i=1,NROW
        DO 160 j=1,NCOL
	   QZEXTRA(j,i,k)=0
	   headdiff=0
	   IF (k.NE.1) THEN
	    DO iz=2,NZONES+1
	     headdiff=headdiff-
     &	     NUS(iz-1)*(ZETA(j,i,k-1,iz)-ZETA(j,i,k-1,iz-1))
          END DO
		  QZEXTRA(j,i,k)=-CV(j,i,k-1)*
     &	    (headdiff+0.5*(ZETA(j,i,k-1,NZONES+1)-ZETA(j,i,k,1))
     &                      *(NUBOT(j,i,k-1)+NUTOP(j,i,k)))
	   END IF
  160 CONTINUE
C
C-------COPY RHS TO RHSFRESH FOR LATER USE IN SWI1BD (SWI 1.0)
      RHSFRESH = RHS
!      DO 175 k=1,NLAY
!       DO 175 i=1,NROW
!        DO 175 j=1,NCOL
!	    RHSFRESH(j,i,k)=RHS(j,i,k)
!  175 CONTINUE

C
C-------ADD QZEXTRA TO RHS
      DO 170 k=1,NLAY
       DO 170 i=1,NROW
        DO 170 j=1,NCOL
	    !EERST=RHS(j,i,k) 
	    RHS(j,i,k)=RHS(j,i,k)+QZEXTRA(j,i,k)
	   IF (k.NE.NLAY) THEN
	    RHS(j,i,k)=RHS(j,i,k)-QZEXTRA(j,i,k+1)
         END IF
		!LATER=RHS(j,i,k) 
	
  170 CONTINUE
C
C------ADD TO RHS DUE TO DENSITY DIFFERENCES
      DO 180 k=1,NLAY
       DO 180 i=1,NROW
        DO 180 j=1,NCOL
	    !EERST=RHS(j,i,k) 

	   DO 180 iz=1,NZONES

C -------LEFT FACE--------------------
		IF (j.NE.1) THEN
	      RHS(j,i,k)=RHS(j,i,k)-DELNUS(iz)*
     &	      SWICUMCR(j-1,i,k,iz)*(ZETA(j-1,i,k,iz)-ZETA(j,i,k,iz))+
     &          EPS(iz)*SWICR(j-1,i,k,iz)*(
     &		  (ZETA(j-1,i,k,iz)-ZETA(j-1,i,k,iz+1))-
     &           (ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)))
	    END IF
    
C -------RIGHT FACE--------------------
		IF (j.NE.NCOL) THEN
	      RHS(j,i,k)=RHS(j,i,k)-DELNUS(iz)*
     &	      SWICUMCR(j,i,k,iz)*(ZETA(j+1,i,k,iz)-ZETA(j,i,k,iz))+
     &          EPS(iz)*SWICR(j,i,k,iz)*(
     &		  (ZETA(j+1,i,k,iz)-ZETA(j+1,i,k,iz+1))-
     &           (ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)))
	    END IF

C -------BACK FACE--------------------	   
		IF (i.NE.1) THEN
	      RHS(j,i,k)=RHS(j,i,k)-DELNUS(iz)*
     &	      SWICUMCC(j,i-1,k,iz)*(ZETA(j,i-1,k,iz)-ZETA(j,i,k,iz))+
     &          EPS(iz)*SWICC(j,i-1,k,iz)*(
     &		  (ZETA(j,i-1,k,iz)-ZETA(j,i-1,k,iz+1))-
     &           (ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)))
	    END IF

C -------FRONT FACE--------------------
		IF (i.NE.NROW) THEN
	      RHS(j,i,k)=RHS(j,i,k)-DELNUS(iz)*
     &	      SWICUMCC(j,i,k,iz)*(ZETA(j,i+1,k,iz)-ZETA(j,i,k,iz))+
     &          EPS(iz)*SWICC(j,i,k,iz)*(
     &		  (ZETA(j,i+1,k,iz)-ZETA(j,i+1,k,iz+1))-
     &           (ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)))

	    END IF
		!LATER=RHS(j,i,k) 
!	   IF (EERST.NE.LATER) THEN
!	    TTTT=3
!         END IF

  180 CONTINUE
C-------RETURN
      RETURN
      END SUBROUTINE GWF2SWI1FM


      SUBROUTINE GWF2SWI1BD(Kkstp,Kkper,Igrid)
C-----VERSION 3 26SEP2010 GWF2SWI1BD
C     ******************************************************************
C     CALCULATE AND SAVE ZETA FOR SWI PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                      LBOTM,BOTM,
     3                      CR,CC,CV,HCOF,RHS,
     4                      DELR,DELC,IBOUND,HNEW,HOLD,
     5                      BUFF,ISSFLG,NSTP
      USE GWFBASMODULE, ONLY: DELT,HDRY
      !USE GWFBCFMODULE, ONLY: SC1
      USE GWFSWIMODULE
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: Kkstp
      INTEGER, INTENT(IN) :: Kkper
      INTEGER, INTENT(IN) :: Igrid
C     + + + LOCAL DEFINITIONS + + +
      INTEGER :: i, j, k
      INTEGER :: iz, iz2, iz3
      INTEGER :: izrev
      INTEGER :: iusezone, iactive
      INTEGER :: icount
      INTEGER :: mxiterde45, niter
      INTEGER :: kkiter
      INTEGER :: iprd4, ifreq, ierr, icnvg
      !INTEGER :: nstp
      REAL :: q, qztop, qzbot
      REAL :: accl, hclose
      REAL :: one, tled, rho
      REAL :: strg
      REAL :: swidelt
      REAL :: zt, zb, zetac, zetaavg
      REAL :: nuontop, nubelbot
	CHARACTER*16 ZETANAME
C
C-------SET POINTERS FOR GRID
      CALL SGWF2SWI1PNT(Igrid)
C
C0-------COPY ZETA TO ZETAOLD
      DO k=1,NLAY
        DO  i=1,NROW
          DO  j=1,NCOL
		    BUFF(j,i,k) = 0.0
	      DO iz = 1, NZONES+1
			    ZETAOLD(j,i,k,iz)=ZETA(j,i,k,iz)
	      END DO
	    END DO
	  END DO
	END DO

C1-------CALCULATE BRHS
      DO 90 k=1,NLAY
       DO 100 i=1,NROW
        DO 100 j=1,NCOL
	   DO 100 iz=1,NZONES
C
C            INITIALIZE BRHS
	     BRHS(j,i,iz)=0.
C
C             CALCULATE BOUNDARY FLUX
		    q = RHSFRESH(j,i,k)-HNEW(j,i,k)*(HCOF(j,i,k)-HCOFOLD(j,i,k)) 
		
            IF(ISSFLG(Kkper).NE.0)  THEN
		      strg=0
!              iusezone=IZONENR(j,i,k)
!	         IF ((IZONENR(j,i,k).LT.0).AND.(q.GT.0)) THEN
!			      iusezone=1
!	         ENDIF
	      ELSE
		      one=1.
		      tled=one/DELT
		      rho=SC1(j,i,k)*tled
		      q = q + rho*HOLD(j,i,k) 
!	         iusezone=IZONENR(j,i,k)
		      strg = rho*HOLD(j,i,k) - rho*HNEW(j,i,k) 
		      q = q - strg
!	         IF ((IZONENR(j,i,k).LT.0).AND.(q.GT.0)) THEN
!		         iusezone=1
!	         ENDIF
	      ENDIF
C             SET ZONE NUMBER FOR BOUNDARY CONDITIONS	
            iusezone=IZONENR(j,i,k)
	      IF ( (IZONENR(j,i,k).LT.0) .AND. (q.GT.0) ) THEN
	        iusezone=1
	      ENDIF

	      IF (IPLPOS(j,i,k,iz).EQ.0) THEN
	        IF ((iz.LE.ABS(iusezone)).AND.(q.NE.0)) THEN
                  BRHS(j,i,iz)=BRHS(j,i,iz) + q
	        ENDIF
	      ELSE 
	        BRHS(j,i,iz)=0
	      ENDIF
C		
	      IF (iz.NE.1) THEN
             IF (IPLPOS(j,i,k,iz).EQ.0) THEN
               BRHS(j,i,iz)= BRHS(j,i,iz) - 
     &	         SSZ(j,i,k)*DELR(j)*DELC(i)*ZETA(j,i,k,iz)/DELT
               CALL SWI1SD(BRHS(1,1,iz),IPLPOS(1,1,k,iz),
     &              SWICUMCR(1,1,k,iz),SWICUMCC(1,1,k,iz),1.,
     &              HNEW(1,1,k),j,i,NCOL,NROW)
               DO iz2=1,iz-1
                 CALL SWI1SR(BRHS(1,1,iz),IPLPOS(1,1,k,iz),
     &                SWICUMCR(1,1,k,iz),SWICUMCC(1,1,k,iz),DELNUS(iz2),
     &                ZETA(1,1,k,iz2),j,i,NCOL,NROW)
               END DO 
             END IF
             CALL SWI1SR(BRHS(1,1,iz),IPLPOS(1,1,k,iz),SWICR(1,1,k,iz),
     &          SWICC(1,1,k,iz),EPS(iz),ZETA(1,1,k,iz+1),j,i,NCOL,NROW)
C
C-----------Zones iz+1 through bottom
             DO iz2 = iz+1, NZONES
               IF (IPLPOS(j,i,k,iz).EQ.0) THEN
                 CALL SWI1SR(BRHS(1,1,iz),IPLPOS(1,1,k,iz),
     &              SWICUMCR(1,1,k,iz2),SWICUMCC(1,1,k,iz2),DELNUS(iz2),
     &              ZETA(1,1,k,iz2),j,i,NCOL,NROW)
               END IF
C----------------LEFT FACE
               IF (j.NE.1) THEN
                 BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &	         SWICR(j-1,i,k,iz2)*((ZETA(j-1,i,k,iz2)-ZETA(j,i,k,iz2))-
     &           (ZETA(j-1,i,k,iz2+1)-ZETA(j,i,k,iz2+1)))
	         END IF
C----------------RIGHT FACE
		       IF (j.NE.NCOL) THEN
                 BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &	         SWICR(j,i,k,iz2)*((ZETA(j+1,i,k,iz2)-ZETA(j,i,k,iz2))-
     &           (ZETA(j+1,i,k,iz2+1)-ZETA(j,i,k,iz2+1)))
	         END IF
C----------------BACK FACE
		       IF (i.NE.1) THEN
                 BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &	         SWICC(j,i-1,k,iz2)*((ZETA(j,i-1,k,iz2)-ZETA(j,i,k,iz2))-
     &           (ZETA(j,i-1,k,iz2+1)-ZETA(j,i,k,iz2+1)))
	         END IF
C----------------FRONT FACE
		       IF (i.NE.NROW) THEN
                 BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &	         SWICC(j,i,k,iz2)*((ZETA(j,i+1,k,iz2)-ZETA(j,i,k,iz2))-
     &	         (ZETA(j,i+1,k,iz2+1)-ZETA(j,i,k,iz2+1)))
	         END IF
	       END DO
C
C2-------------SPECIFY BOUNDARY CONDITION OF FLOW AT TIPS AND TOES
  	       IF (IPLPOS(j,i,k,iz).EQ.0) THEN
C----------------LEFT FACE
               IF (j.NE.1) THEN
		         IF (IPLPOS(j-1,i,k,iz).NE.0) THEN
                   CALL SWI1SQR (q,HNEW,ZETA,NZONES,DELNUS,EPS,
     &                  SWICUMCR,SWICR,j,i,k,iz,NCOL,NROW,NLAY)
	             BRHS(j,i,iz)=BRHS(j,i,iz)-q
     	           END IF
	         END IF
C----------------RIGHT FACE
		       IF (j.NE.NCOL) THEN
		         IF (IPLPOS(j+1,i,k,iz).NE.0) THEN
                   CALL SWI1SQR (q,HNEW,ZETA,NZONES,DELNUS,EPS,
     &                  SWICUMCR,SWICR,j+1,i,k,iz,NCOL,NROW,NLAY)
	             BRHS(j,i,iz)=BRHS(j,i,iz)+q
		         END IF
	         END IF
C----------------BACK FACE
    		       IF (i.NE.1) THEN
		         IF (IPLPOS(j,i-1,k,iz).NE.0) THEN
                   CALL SWI1SQC (q,HNEW,ZETA,NZONES,DELNUS,EPS,
     &                  SWICUMCC,SWICC,j,i,k,iz,NCOL,NROW,NLAY)
	             BRHS(j,i,iz)=BRHS(j,i,iz)-q
	           END IF
	         END IF
C----------------FRONT FACE
		       IF (i.NE.NROW) THEN
		         IF (IPLPOS(j,i+1,k,iz).NE.0) THEN
                   CALL SWI1SQC (q,HNEW,ZETA,NZONES,DELNUS,EPS,
     &                  SWICUMCC,SWICC,j,i+1,k,iz,NCOL,NROW,NLAY)
	             BRHS(j,i,iz)=BRHS(j,i,iz)+q
		         END IF
	         END IF
C----------------UPPER FACE
		       IF (k.NE.1) THEN
                 qztop=CV(j,i,k-1)*(HNEW(j,i,k)-HNEW(j,i,k-1))+
     &                 QZEXTRA(j,i,k)
  			       nuontop=NUBOT(j,i,k-1)
			       IF ((qztop.LT.0).AND.(nuontop.GE.NUS(iz)).AND.
     &		           (NUBOT(j,i,k).GE.nuontop)) THEN
		           BRHS(j,i,iz)=BRHS(j,i,iz)+qztop
			       ENDIF
	         END IF
C----------------LOWER FACE
		       IF (k.NE.NLAY) THEN
                 qzbot=CV(j,i,k)*(HNEW(j,i,k)-HNEW(j,i,k+1))-
     &                 QZEXTRA(j,i,k+1)
	           nubelbot=NUTOP(j,i,k+1)	
			       IF ((qzbot.LT.0).AND.(nubelbot.LT.NUS(iz)).AND.
     &		           (NUTOP(j,i,k).LE.nubelbot)) THEN
                   BRHS(j,i,iz)=BRHS(j,i,iz)+0
			       ELSE
                   BRHS(j,i,iz)=BRHS(j,i,iz)+qzbot
                 END IF
	         END IF
   	       END IF
	     END IF

C	BUFF(j,i,iz)=BRHS(j,i,iz)

  100 CONTINUE
C
C3----- CALL DE45 TO SOLVE EQUATIONS FOR 1 LAYER
	!iactive=0 !COMMENTED OUT MF2005 VERSION 3.0
	DO 115 iz=2,NZONES
   	 iactive=0  !MOVED FROM ABOVE FOR MF2005 VERSION 3.0
       DO 110 i=1,NROW
        DO 110 j=1,NCOL
	   IF (j.NE.NCOL) THEN
         IF ((IPLPOS(j,i,k,iz).EQ.0).AND.(IPLPOS(j+1,i,k,iz).EQ.0)) THEN
	      SWISOLVCR(j,i)=DELNUS(iz)*SWICUMCR(j,i,k,iz)
         ELSE  
	      SWISOLVCR(j,i)=0.
	   END IF
	   END IF
         SWISOLVCR(j,i)=SWISOLVCR(j,i)-EPS(iz)*SWICR(j,i,k,iz)

	   IF (i.NE.NROW) THEN
         IF ((IPLPOS(j,i,k,iz).EQ.0).AND.(IPLPOS(j,i+1,k,iz).EQ.0)) THEN
	      SWISOLVCC(j,i)=DELNUS(iz)*SWICUMCC(j,i,k,iz)
         ELSE  
	      SWISOLVCC(j,i)=0.
	   END IF
	   END IF
         SWISOLVCC(j,i)=SWISOLVCC(j,i)-EPS(iz)*SWICC(j,i,k,iz)

	   SWISOLVCV(j,i)=0
	   SWIHCOF(j,i)=-SSZ(j,i,k)*DELR(j)*DELC(i)/DELT

	   DUM(j,i)=ZETA(j,i,k,iz)

         IF (IPLPOS(j,i,k,iz).EQ.0) THEN
	    IBO(j,i)=1
	    iactive = iactive + 1
	   ELSE
          IBO(j,i)=0
	   END IF
  110 CONTINUE
C
C-------SOLVE FOR ZETA SURFACE
	IF (USEDE45.EQ.1) THEN
C5--------INITIALIZE PARAMETERS FOR DIRECT SOLVING
!        ITMX=1
!        MXITER=ITMX
!	  niter=1
!	  accl=1.0
!	  hclose=0.00001
!	  iprd4=0
!	  ifreq=3
	  kkiter=1
	  swidelt = SWIDE4%DELTL
      	ierr=0

!      CALL ULASAV(BUFF(1,1,1),'BRHS1           ',
!     & kkstp,Kkper,0.0,0.0,NCOL,
!     &                NROW,k,333)
!      CALL ULASAV(BUFF(1,1,2),'BRHS2            ',
!     & kkstp,Kkper,0.0,0.0,NCOL,
!     &                NROW,k,444)


C6-------CALL DE45AP TO SOLVE FOR ZETA
!      SUBROUTINE DE47AP(HNEW,IBOUND,AU,AL,IUPPNT,IEQPNT,D4B,MXUP,MXLOW,
!     1  MXEQ,MXBW,CR,CC,CV,HCOF,RHS,ACCLDE4,KITER,ITMX,MXITER,NITERDE4,
!     2  HCLOSEDE4,IPRD4,ICNVG,NCOL,NROW,NLAY,IOUT,LRCHDE4,HDCGDE4,
!     3  IFREQ,KSTP,KPER,DELT,NSTP,ID4DIR,ID4DIM,MUTD4,DELTL,NBWL,NUPL,
!     4  NLOWL,NLOW,NEQ,NUP,NBW,IERR)

      CALL DE47AP(DUM,IBO,SWIDE4%AU,SWIDE4%AL,SWIDE4%IUPPNT,
     2  SWIDE4%IEQPNT,SWIDE4%D4B,SWIDE4%MXUP,SWIDE4%MXLOW,
     3  SWIDE4%MXEQ,SWIDE4%MXBW,SWISOLVCR,SWISOLVCC,SWISOLVCV,SWIHCOF,
     4  BRHS(1,1,iz),SWIDE4%ACCLDE4,kkiter,SWIDE4%ITMX,SWIDE4%MXITER,
     5  SWIDE4%NITERDE4,SWIDE4%HCLOSEDE4,SWIDE4%IPRD4,icnvg,
     6  NCOL,NROW,NLAYSWI,IOUT,SWIDE4%LRCHDE4,SWIDE4%HDCGDE4,
     7  SWIDE4%IFREQ,Kkstp,Kkper,swidelt,NSTP(Kkper),SWIDE4%ID4DIR,
     8  SWIDE4%ID4DIM,SWIDE4%MUTD4,SWIDE4%DELTL,SWIDE4%NBWL,SWIDE4%NUPL,
     9  SWIDE4%NLOWL,SWIDE4%NLOW,SWIDE4%NEQ,SWIDE4%NUP,SWIDE4%NBW,ierr)

!        CALL DE47AP(DUM,IBO,AU,AL,IUPP,IEQP,SWID4B,MXUP,MXLOW,
!     &              MXEQ,MXBW,SWISOLVCR,SWISOLVCC,SWISOLVCV,SWIHCOF,
!     &              BRHS(1,1,iz),accl,kkiter,ITMX,MXITER,niter,
!     &              hclose,iprd4,icnvg,NCOL,NROW,NLAYSWI,IOUT,
!     &              SWILRCH,SWIHDCG,
!     &              ifreq,Kkstp,Kkper,swidelt,NSTP(Kkper),
!     &              ID4DIR,ID4DIM,MUTD4,DELTL,NBWL,NUPL,
!     &              NLOWL,NLOW,NEQ,NUP,NBW,ierr)

	  IF (ierr.EQ.0) THEN
	    WRITE(IOUT,'(A15,I3,A10,I3,A6,I3,A6,I3,A15,I5,A20,I8)') 
     &	    ' STRESS PERIOD=',Kkper,
     &	    ' TIMESTEP=',kkstp,
     &	    ' LAYER=',k,
     &	    ' ZONE=',iz,
     &	    ' ACTIVE CELLS=',iactive,
     &	    ' NUMBER OF EQUATIONS=',SWIDE4%NEQ  !ierr	   
	  END IF
C-------PCG
	ELSE
!	  NODES=NCOL*NROW
!        MXITER=ITMX
!	  niter=1
!	  accl=1.0
!	  hclose=0.00001
!	  iprd4=0
!	  ifreq=3
!	  swidelt=1.
	  kkiter=1
      	ierr = 0

!      SUBROUTINE PCG7AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,V,SS,P,CD,HCHG,
!     &                  LHCH,RCHG,LRCH,KITER,NITER,HCLOSE,RCLOSE,
!     &                  ICNVG,KSTP,KPER,IPRPCG,MXITER,ITER1,NPCOND,
!     &                  NBPOL,NSTP,NCOL,NROW,NLAY,NODES,RELAX,IOUT,
!     &                  MUTPCG,IT1,DAMPSS,RES,HCSV,IERR,HPCG,DAMPTR,
!     &                  ISS,HDRY)

      CALL PCG7AP(DUM,IBO,SWISOLVCR,SWISOLVCC,SWISOLVCV,SWIHCOF,
     2  BRHS(1,1,iz),SWIPCG%VPCG,SWIPCG%SS,SWIPCG%P,SWIPCG%CD,
     3  SWIPCG%HCHG,SWIPCG%LHCH,SWIPCG%RCHG,SWIPCG%LRCHPCG,kkiter,
     4  SWIPCG%NITER,SWIPCG%HCLOSEPCG,SWIPCG%RCLOSEPCG,icnvg,
     5  Kkstp,Kkper,SWIPCG%IPRPCG,SWIPCG%MXITER,SWIPCG%ITER1,
     6  SWIPCG%NPCOND,SWIPCG%NBPOL,NSTP(Kkper),NCOL,NROW,NLAYSWI,
     7  SWIPCG%NODES,SWIPCG%RELAXPCG,IOUT,
     7  SWIPCG%MUTPCG,SWIPCG%IT1,SWIPCG%DAMPPCG,BUFF,SWIPCG%HCSV,
     8  ierr,SWIPCG%HPCG,SWIPCG%DAMPPCGT,ISSFLG(Kkper),HDRY)


!            CALL  PCG2AP(DUM,IBO,SWISOLVCR,SWISOLVCC,
!     &              SWISOLVCV,SWIHCOF,BRHS(1,1,iz),V,
!     &               SS,P,CD,HCHG,
!     &               LHCH,RCHG,LRCH,kkiter,
!     &              niter,hclose,RCLOSE,icnvg,kkstp,Kkper,
!     &               IPRPCG,MXITER,ITER1,NPCOND,NBPOL,
!     &                nstp,NCOL,NROW,1,NODES,RELAX,
!     &               IOUT,MUTPCG,IT1,DAMP,BUFF,
!     &               HCSV,ierr,ierru,HPCG)


	END IF



!CDEBUG
!	DO  i=1,NROW
!        DO  j=1,NCOL
!	   BUFF(j,i,1)=IBO(j,i)
!	   ENDDO
!	ENDDO	
!
!      CALL ULASAV(BUFF(1,1,1),'IBO             ',
!     & kkstp,Kkper,0.0,0.0,NCOL,
!     &                NROW,k,555)
!
!CDEBUG



C
C7-------COPY REAL ARRAY DUM INTO DOUBLE ARRAY ZETA
       ZETA(:,:,k,iz) = DUM(:,:)
!       DO 95 i=1,NROW
!        DO 95 j=1,NCOL
!	   ZETA(j,i,k,iz)=DUM(j,i)
!   95 CONTINUE
  115 CONTINUE
   90 CONTINUE
C
C8-----MOVE TIPS AND TOES
C8-----NOTE: MAY DIFFER FROM MATLAB PROTOTYPE BECAUSE ZETA IS 
C8-----ALTERED DURING THE SEARCH PROCESS
      DO 120 k=1,NLAY
	 DO 120 iz=2,NZONES
        DO 120 i=1,NROW
         DO 120 j=1,NCOL

	    zt=ZETA(j,i,k,1)
	    zb=ZETA(j,i,k,NZONES+1)
	    zetac=ZETA(j,i,k,iz)



	   IF (IPLPOS(j,i,k,iz).EQ.0) THEN
          IF ((j.NE.1).AND.(j.NE.NCOL)) THEN
C -------LEFT FACE--------------------
		 IF (IPLPOS(j-1,i,k,iz).EQ.1) THEN
            IF ((zt-zetac).GT.(TIPSLOPE*0.5*(DELR(j)+DELR(j-1)))) THEN
             ZETA(j,i,k,iz)=zetac+DELZETA*DELR(j-1)/DELR(j)
             ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,1)-DELZETA
	      ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
	       IF (IPLPOS(j+1,i,k,iz).EQ.0) THEN
              ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,iz)-(zt-zetac)
              ZETA(j,i,k,iz)=zt
	       END IF
            END IF
           ELSEIF (IPLPOS(j-1,i,k,iz).EQ.2) THEN
            IF (zetac.GT.(zb+TOESLOPE*0.5*(DELR(j)+DELR(j-1)))) THEN
             ZETA(j,i,k,iz)=zetac-DELZETA*DELR(j-1)/DELR(j)
             ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,NZONES+1)+DELZETA
	      ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
	       IF (IPLPOS(j+1,i,k,iz).EQ.0) THEN
              ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,iz)+(zetac-zb)
              ZETA(j,i,k,iz)=zb
	       END IF
            END IF
     	     END IF
C -------RIGHT FACE--------------------
		 IF (IPLPOS(j+1,i,k,iz).EQ.1) THEN
            IF ((zt-zetac).GT.(TIPSLOPE*0.5*(DELR(j)+DELR(j+1)))) THEN
             ZETA(j,i,k,iz)=zetac+DELZETA*DELR(j+1)/DELR(j)
             ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,1)-DELZETA
	      ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
	       IF (IPLPOS(j-1,i,k,iz).EQ.0) THEN
              ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,iz)-(zt-zetac)
              ZETA(j,i,k,iz)=zt
	       END IF
            END IF
           ELSEIF (IPLPOS(j+1,i,k,iz).EQ.2) THEN
            IF (zetac.GT.(zb+TOESLOPE*0.5*(DELR(j)+DELR(j+1)))) THEN
             ZETA(j,i,k,iz)=zetac-DELZETA*DELR(j+1)/DELR(j)
             ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,NZONES+1)+DELZETA
	      ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
		   IF (IPLPOS(j-1,i,k,iz).EQ.0) THEN
              ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,iz)+(zetac-zb)
              ZETA(j,i,k,iz)=zb
	       END IF
            END IF
     	     END IF

	    END IF

          IF ((i.NE.1).AND.(i.NE.NROW)) THEN
C -------BACK FACE--------------------
		 IF (IPLPOS(j,i-1,k,iz).EQ.1) THEN
            IF ((zt-zetac).GT.(TIPSLOPE*0.5*(DELC(i)+DELC(i-1)))) THEN
             ZETA(j,i,k,iz)=zetac+DELZETA*DELC(i-1)/DELC(i)
             ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,1)-DELZETA
	      ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
		   IF (IPLPOS(j,i+1,k,iz).EQ.0) THEN
              ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,iz)-(zt-zetac)
              ZETA(j,i,k,iz)=zt
	       END IF
            END IF
           ELSEIF (IPLPOS(j,i-1,k,iz).EQ.2) THEN
            IF (zetac.GT.(zb+TOESLOPE*0.5*(DELC(i)+DELC(i-1)))) THEN
             ZETA(j,i,k,iz)=zetac-DELZETA*DELC(i-1)/DELC(i)
             ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,NZONES+1)+DELZETA
	      ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
		   IF (IPLPOS(j,i+1,k,iz).EQ.0) THEN
              ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,iz)+(zetac-zb)
              ZETA(j,i,k,iz)=zb
	       END IF
            END IF
     	     END IF
C -------FRONT FACE--------------------
		 IF (IPLPOS(j,i+1,k,iz).EQ.1) THEN
            IF ((zt-zetac).GT.(TIPSLOPE*0.5*(DELC(i)+DELC(i+1)))) THEN
             ZETA(j,i,k,iz)=zetac+DELZETA*DELC(i+1)/DELC(i)
             ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,1)-DELZETA
	      ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
	       IF (IPLPOS(j,i-1,k,iz).EQ.0) THEN
              ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,iz)-(zt-zetac)
              ZETA(j,i,k,iz)=zt
	       END IF
            END IF
           ELSEIF (IPLPOS(j,i+1,k,iz).EQ.2) THEN
            IF (zetac.GT.(zb+TOESLOPE*0.5*(DELC(i)+DELC(i+1)))) THEN
             ZETA(j,i,k,iz)=zetac-DELZETA*DELC(i+1)/DELC(i)
             ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,NZONES+1)+DELZETA
	      ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
	       IF (IPLPOS(j,i-1,k,iz).EQ.0) THEN
              ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,iz)+(zetac-zb)
              ZETA(j,i,k,iz)=zb
	       END IF
            END IF
     	     END IF

	    END IF
      
         END IF


  120 CONTINUE
C
C9------CHECK WHETHER ANYWHERE THE THICKNESS GETS TOO THIN
      DO 130 k=1,NLAY
       DO 130 i=1,NROW
        DO 130 j=1,NCOL
	   DO 130 iz=2,NZONES
	    IF (IPLPOS(j,i,k,iz).EQ.0) THEN
           IF ((ZETA(j,i,k,iz)-ZETA(j,i,k,NZONES+1)).LT.ZETAMIN) THEN
              ZETA(j,i,k,iz)=ZETA(j,i,k,NZONES+1)
           END IF
           IF ((ZETA(j,i,k,1)-ZETA(j,i,k,iz)).LT.ZETAMIN) THEN
              ZETA(j,i,k,iz)=ZETA(j,i,k,1)
           END IF
	    END IF
  130 CONTINUE 	
C
C10------FIND IF ANYWHERE THE PLANES ARE CROSSING
      DO 150 k=1,NLAY
       DO 150 i=1,NROW
        DO 150 j=1,NCOL
	   DO 150 iz=2,NZONES-1
          IF ((ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)).LT.0.001) THEN
             zetaavg=0.5*(ZETA(j,i,k,iz)+ZETA(j,i,k,iz+1))
             ZETA(j,i,k,iz)=zetaavg
             ZETA(j,i,k,iz+1)=zetaavg
		   DO iz2=2,iz-1
              izrev=iz+1-iz2
              IF ((ZETA(j,i,k,izrev)-ZETA(j,i,k,iz+1)).LT.0.001) THEN
               zetaavg=0.
	         icount=0
               DO iz3=izrev,iz+1
                icount=icount+1
                zetaavg=zetaavg + ZETA(j,i,k,iz3)
	         END DO
               IF (icount.NE.0) THEN
                zetaavg=zetaavg/icount
                DO iz3=izrev,iz+1
                 ZETA(j,i,k,iz3)=zetaavg
                END DO
	         END IF
	        END IF
	       END DO
          END IF
  150 CONTINUE
C            
C11------WRITE ZETA TO UNFORMATTED FILE
	IF (ISWIZT.NE.0) THEN
	 IF (NPRS.GE.0) THEN
	  DO iz=2,NZONES
         WRITE(ZETANAME,151) iz-1
 151     FORMAT('    ZETAPLANE ',I2)
	   IF(MOD(kkstp,NPRS).EQ.0) THEN
C	   WRITE (ISWIZT) ZETA
         CALL UBUDSV(kkstp,Kkper,ZETANAME,
     &                ISWIZT,ZETA(1,1,1,iz),NCOL,NROW,NLAY,IOUT) 

	         CALL UBUDSV(kkstp,Kkper,'QZE             ',
     &                222,QZEXTRA,NCOL,NROW,NLAY,IOUT) 


         ELSEIF ((NPRS.EQ.0).AND.(kkstp.EQ.NSTP(Kkper))) THEN
C	    WRITE (ISWIZT) ZETA
         CALL UBUDSV(kkstp,Kkper,ZETANAME,
     &                ISWIZT,ZETA(1,1,1,iz),NCOL,NROW,NLAY,IOUT) 
	         CALL UBUDSV(kkstp,Kkper,'QZE             ',
     &                222,QZEXTRA,NCOL,NROW,NLAY,IOUT) 
         END IF
	  END DO
	 END IF
	END IF
C      
C12-----RETURN
      RETURN
      END SUBROUTINE GWF2SWI1BD


      SUBROUTINE SWI1SQR (Q,HNEW,ZETA,NZONES,DELNUS,EPS,SWICUMCR,SWICR,
     &                                        J,I,K,IZ,NCOL,NROW,NLAY)
     

C     ******************************************************************
C     SWI
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
	DIMENSION ZETA(NCOL,NROW,NLAY,NZONES+1)
	DOUBLE PRECISION HNEW
      DIMENSION HNEW(NCOL,NROW,NLAY)
	DIMENSION SWICR(NCOL,NROW,NLAY,NZONES)
C	DIMENSION SWICC(NCOL,NROW,NLAY,NZONES)
      DIMENSION SWICUMCR(NCOL,NROW,NLAY,NZONES)
C	DIMENSION SWICUMCC(NCOL,NROW,NLAY,NZONES)
	REAL NUS(NZONES)
      REAL EPS(NZONES)
      REAL DELNUS(NZONES)
C
C1-----CALCULATE Q IN ROW DIRECTION      
      Q=SWICUMCR(J-1,I,K,IZ)*(HNEW(J-1,I,K)-HNEW(J,I,K))

	DO IZ2=1,IZ
       Q=Q+SWICUMCR(J-1,I,K,IZ)*DELNUS(IZ)*
     &                 (ZETA(J-1,I,K,IZ2)-ZETA(J,I,K,IZ2))
	END DO

      Q=Q-SWICR(J-1,I,K,IZ)*EPS(IZ)*
     &                 (ZETA(J-1,I,K,IZ)-ZETA(J,I,K,IZ))
      Q=Q+SWICR(J-1,I,K,IZ)*EPS(IZ)*
     &                 (ZETA(J-1,I,K,IZ+1)-ZETA(J,I,K,IZ+1))


      DO IZ2=IZ+1,NZONES
       Q=Q+SWICUMCR(J-1,I,K,IZ2)*DELNUS(IZ2)*
     &                 (ZETA(J-1,I,K,IZ2)-ZETA(J,I,K,IZ2))
       Q=Q-SWICR(J-1,I,K,IZ2)*EPS(IZ2)*
     &                 (ZETA(J-1,I,K,IZ2)-ZETA(J,I,K,IZ2))
       Q=Q+SWICR(J-1,I,K,IZ2)*EPS(IZ2)*
     &                 (ZETA(J-1,I,K,IZ2+1)-ZETA(J,I,K,IZ2+1))
	END DO
C     ------------------------------------------------------------------
C
C2------RETURN
      RETURN
      END SUBROUTINE SWI1SQR



      SUBROUTINE SWI1SQC (Q,HNEW,ZETA,NZONES,DELNUS,EPS,SWICUMCC,SWICC,
     &                                        J,I,K,IZ,NCOL,NROW,NLAY)
     

C     ******************************************************************
C     SWI
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
	DIMENSION ZETA(NCOL,NROW,NLAY,NZONES+1)
	DOUBLE PRECISION HNEW 
      DIMENSION HNEW(NCOL,NROW,NLAY)
C	DIMENSION SWICR(NCOL,NROW,NLAY,NZONES)
	DIMENSION SWICC(NCOL,NROW,NLAY,NZONES)
C     DIMENSION SWICUMCR(NCOL,NROW,NLAY,NZONES)
	DIMENSION SWICUMCC(NCOL,NROW,NLAY,NZONES)
	REAL NUS(NZONES),EPS(NZONES),DELNUS(NZONES)
C
C1-----CALCULATE Q IN COLUMN DIRECTION      
      Q=SWICUMCC(J,I-1,K,IZ)*(HNEW(J,I-1,K)-HNEW(J,I,K))

	DO IZ2=1,IZ
       Q=Q+SWICUMCC(J,I-1,K,IZ)*DELNUS(IZ)*
     &                 (ZETA(J,I-1,K,IZ2)-ZETA(J,I,K,IZ2))
	END DO

      Q=Q-SWICC(J,I-1,K,IZ)*EPS(IZ)*
     &                 (ZETA(J,I-1,K,IZ)-ZETA(J,I,K,IZ))
      Q=Q+SWICC(J,I-1,K,IZ)*EPS(IZ)*
     &                 (ZETA(J,I-1,K,IZ+1)-ZETA(J,I,K,IZ+1))


      DO IZ2=IZ+1,NZONES
       Q=Q+SWICUMCC(J,I-1,K,IZ2)*DELNUS(IZ2)*
     &                 (ZETA(J,I-1,K,IZ2)-ZETA(J,I,K,IZ2))
       Q=Q-SWICC(J,I-1,K,IZ2)*EPS(IZ2)*
     &                 (ZETA(J,I-1,K,IZ2)-ZETA(J,I,K,IZ2))
       Q=Q+SWICC(J,I-1,K,IZ2)*EPS(IZ2)*
     &                 (ZETA(J,I-1,K,IZ2+1)-ZETA(J,I,K,IZ2+1))
	END DO
C     ------------------------------------------------------------------
C
C2------RETURN
      RETURN
      END SUBROUTINE SWI1SQC



      SUBROUTINE SWI1SD(B,IPOS,CRLAY,CCLAY,FAC,VAR,J,I,NCOL,NROW)

C     ******************************************************************
C     SWI
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
	DIMENSION B(NCOL,NROW),CRLAY(NCOL,NROW),CCLAY(NCOL,NROW)
	DIMENSION IPOS(NCOL,NROW)
	DOUBLE PRECISION VAR(NCOL,NROW)
C
C1-----CALCULATE FLUX IN COLUMN DIRECTION      
C -------LEFT FACE--------------------
          IF (J.NE.1) THEN
		 IF (IPOS(J-1,I).EQ.0) THEN
            B(J,I) = B(J,I)- FAC*CRLAY(J-1,I)*(VAR(J-1,I)-VAR(J,I))            
	     END IF
	    END IF
    
C -------RIGHT FACE--------------------
		IF (J.NE.NCOL) THEN
		 IF (IPOS(J+1,I).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CRLAY(J,I)*(VAR(J+1,I)-VAR(J,I))
	     END IF
	    END IF

C
C2-----CALCULATE FLUX IN ROW DIRECTION      

C -------BACK FACE--------------------	   
		IF (I.NE.1) THEN
		 IF (IPOS(J,I-1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I-1)*(VAR(J,I-1)-VAR(J,I))
	     END IF
	    END IF

C -------FRONT FACE--------------------
		IF (I.NE.NROW) THEN
		 IF (IPOS(J,I+1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I)*(VAR(J,I+1)-VAR(J,I))
	     END IF
	    END IF

C3------RETURN
      RETURN
      END SUBROUTINE SWI1SD

      
      SUBROUTINE SWI1SR(B,IPOS,CRLAY,CCLAY,FAC,VAR,J,I,NCOL,NROW)
     

C     ******************************************************************
C     SWI
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
	DIMENSION B(NCOL,NROW),CRLAY(NCOL,NROW),CCLAY(NCOL,NROW)
	DIMENSION IPOS(NCOL,NROW)
	REAL VAR(NCOL,NROW)
C
C1-----CALCULATE FLUX IN COLUMN DIRECTION      
C -------LEFT FACE--------------------
          IF (J.NE.1) THEN
		 IF (IPOS(J-1,I).EQ.0) THEN
            B(J,I) = B(J,I)- FAC*CRLAY(J-1,I)*(VAR(J-1,I)-VAR(J,I))            
	     END IF
	    END IF
    
C -------RIGHT FACE--------------------
		IF (J.NE.NCOL) THEN
		 IF (IPOS(J+1,I).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CRLAY(J,I)*(VAR(J+1,I)-VAR(J,I))
	     END IF
	    END IF

C
C2-----CALCULATE FLUX IN ROW DIRECTION      

C -------BACK FACE--------------------	   
		IF (I.NE.1) THEN
		 IF (IPOS(J,I-1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I-1)*(VAR(J,I-1)-VAR(J,I))
	     END IF
	    END IF

C -------FRONT FACE--------------------
		IF (I.NE.NROW) THEN
		 IF (IPOS(J,I+1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I)*(VAR(J,I+1)-VAR(J,I))
	     END IF
	    END IF

C3------RETURN
      RETURN
      END SUBROUTINE SWI1SR




      SUBROUTINE GWF2SWI1DA(Igrid)
C  Deallocate SWR data for a grid
      USE GWFSWIMODULE
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: Igrid
C
      DEALLOCATE(GWFSWIDAT(Igrid)%NPLN)
      DEALLOCATE(GWFSWIDAT(Igrid)%ISTRAT)
      DEALLOCATE(GWFSWIDAT(Igrid)%ISWIZT)
      DEALLOCATE(GWFSWIDAT(Igrid)%NPRS)
      DEALLOCATE(GWFSWIDAT(Igrid)%NZONES)
      DEALLOCATE(GWFSWIDAT(Igrid)%NLAYSWI)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NODES)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NHALFU)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NHALFL)
      DEALLOCATE(GWFSWIDAT(Igrid)%USEDE45)
!      DEALLOCATE(GWFSWIDAT(Igrid)%ID4DIM)
!      DEALLOCATE(GWFSWIDAT(Igrid)%ID4DIR)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NBWGRD)
!      DEALLOCATE(GWFSWIDAT(Igrid)%ITMX)
!      DEALLOCATE(GWFSWIDAT(Igrid)%MXITER)
!      DEALLOCATE(GWFSWIDAT(Igrid)%MXUP)
!      DEALLOCATE(GWFSWIDAT(Igrid)%MXLOW)
!      DEALLOCATE(GWFSWIDAT(Igrid)%MXBW)
!      DEALLOCATE(GWFSWIDAT(Igrid)%MXEQ)
!      DEALLOCATE(GWFSWIDAT(Igrid)%MUTD4)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NBWL)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NUPL)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NLOWL)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NLOW)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NEQ)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NUP)
!      DEALLOCATE(GWFSWIDAT(Igrid)%NBW)
!      DEALLOCATE(GWFSWIDAT(Igrid)%DELTL)
      DEALLOCATE(GWFSWIDAT(Igrid)%TOESLOPE)
      DEALLOCATE(GWFSWIDAT(Igrid)%TIPSLOPE)
      DEALLOCATE(GWFSWIDAT(Igrid)%ZETAMIN)
      DEALLOCATE(GWFSWIDAT(Igrid)%DELZETA)
!      DEALLOCATE(GWFSWIDAT(Igrid)%AU)
!      DEALLOCATE(GWFSWIDAT(Igrid)%IUPP)
!      DEALLOCATE(GWFSWIDAT(Igrid)%AL)
!      DEALLOCATE(GWFSWIDAT(Igrid)%IEQP)
!      DEALLOCATE(GWFSWIDAT(Igrid)%SWID4B)
!      DEALLOCATE(GWFSWIDAT(Igrid)%SWILRCH)
!      DEALLOCATE(GWFSWIDAT(Igrid)%SWIHDCG)
      DEALLOCATE(GWFSWIDAT(Igrid)%IBO)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWIHCOF)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWISOLVCR)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWISOLVCC)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWISOLVCV)
      DEALLOCATE(GWFSWIDAT(Igrid)%ZETA)
      DEALLOCATE(GWFSWIDAT(Igrid)%ZETAOLD)
      DEALLOCATE(GWFSWIDAT(Igrid)%SSZ)
      DEALLOCATE(GWFSWIDAT(Igrid)%EPS)
      DEALLOCATE(GWFSWIDAT(Igrid)%NUS)
      DEALLOCATE(GWFSWIDAT(Igrid)%DELNUS)
      DEALLOCATE(GWFSWIDAT(Igrid)%NUPLANE)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWICR)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWICC)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWICUMCR)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWICUMCC)
      DEALLOCATE(GWFSWIDAT(Igrid)%NUTOP)
      DEALLOCATE(GWFSWIDAT(Igrid)%NUBOT)
      DEALLOCATE(GWFSWIDAT(Igrid)%BRHS)
      DEALLOCATE(GWFSWIDAT(Igrid)%QZEXTRA)
      DEALLOCATE(GWFSWIDAT(Igrid)%DUM)
      DEALLOCATE(GWFSWIDAT(Igrid)%RHSFRESH)
      DEALLOCATE(GWFSWIDAT(Igrid)%HCOFOLD)
      DEALLOCATE(GWFSWIDAT(Igrid)%IPLPOS)
      DEALLOCATE(GWFSWIDAT(Igrid)%IZONENR)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWIDE4)
      DEALLOCATE(GWFSWIDAT(Igrid)%SWIPCG)
C-------RETURN
      RETURN
      END SUBROUTINE GWF2SWI1DA

      SUBROUTINE SGWF2SWI1PNT(Igrid)
C  Set pointers to SWI data for a grid
      USE GWFSWIMODULE
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: Igrid
C
        NPLN=>GWFSWIDAT(Igrid)%NPLN
        ISTRAT=>GWFSWIDAT(Igrid)%ISTRAT
        ISWIZT=>GWFSWIDAT(Igrid)%ISWIZT
        NPRS=>GWFSWIDAT(Igrid)%NPRS
        NZONES=>GWFSWIDAT(Igrid)%NZONES
        NLAYSWI=>GWFSWIDAT(Igrid)%NLAYSWI
!        NODES=>GWFSWIDAT(Igrid)%NODES
!        NHALFU=>GWFSWIDAT(Igrid)%NHALFU
!        NHALFU=>GWFSWIDAT(Igrid)%NHALFL
        USEDE45=>GWFSWIDAT(Igrid)%USEDE45
!        ID4DIM=>GWFSWIDAT(Igrid)%ID4DIM
!        ID4DIR=>GWFSWIDAT(Igrid)%ID4DIR
!        NBWGRD=>GWFSWIDAT(Igrid)%NBWGRD
!        ITMX=>GWFSWIDAT(Igrid)%ITMX
!        MXITER=>GWFSWIDAT(Igrid)%MXITER
!        MXUP=>GWFSWIDAT(Igrid)%MXUP
!        MXLOW=>GWFSWIDAT(Igrid)%MXLOW
!        MXBW=>GWFSWIDAT(Igrid)%MXBW
!        MXEQ=>GWFSWIDAT(Igrid)%MXEQ
!        MUTD4=>GWFSWIDAT(Igrid)%MUTD4
!        NBWL=>GWFSWIDAT(Igrid)%NBWL
!        NUPL=>GWFSWIDAT(Igrid)%NUPL
!        NLOWL=>GWFSWIDAT(Igrid)%NLOWL
!        NLOW=>GWFSWIDAT(Igrid)%NLOW
!        NEQ=>GWFSWIDAT(Igrid)%NEQ
!        NUP=>GWFSWIDAT(Igrid)%NUP
!        NBW=>GWFSWIDAT(Igrid)%NBW
!        DELTL=>GWFSWIDAT(Igrid)%DELTL
        TOESLOPE=>GWFSWIDAT(Igrid)%TOESLOPE
        TIPSLOPE=>GWFSWIDAT(Igrid)%TIPSLOPE
        ZETAMIN=>GWFSWIDAT(Igrid)%ZETAMIN
        DELZETA=>GWFSWIDAT(Igrid)%DELZETA
!        AU=>GWFSWIDAT(Igrid)%AU
!        IUPP=>GWFSWIDAT(Igrid)%IUPP
!        AL=>GWFSWIDAT(Igrid)%AL
!        IEQP=>GWFSWIDAT(Igrid)%IEQP
!        SWID4B=>GWFSWIDAT(Igrid)%SWID4B
!        SWILRCH=>GWFSWIDAT(Igrid)%SWILRCH
!        SWIHDCG=>GWFSWIDAT(Igrid)%SWIHDCG
        IBO=>GWFSWIDAT(Igrid)%IBO
        SWIHCOF=>GWFSWIDAT(Igrid)%SWIHCOF
        SWISOLVCR=>GWFSWIDAT(Igrid)%SWISOLVCR
        SWISOLVCC=>GWFSWIDAT(Igrid)%SWISOLVCC
        SWISOLVCV=>GWFSWIDAT(Igrid)%SWISOLVCV
        ZETA=>GWFSWIDAT(Igrid)%ZETA
        ZETAOLD=>GWFSWIDAT(Igrid)%ZETAOLD
        SSZ=>GWFSWIDAT(Igrid)%SSZ
        EPS=>GWFSWIDAT(Igrid)%EPS
        NUS=>GWFSWIDAT(Igrid)%NUS
        DELNUS=>GWFSWIDAT(Igrid)%DELNUS
        NUPLANE=>GWFSWIDAT(Igrid)%NUPLANE
        SWICR=>GWFSWIDAT(Igrid)%SWICR
        SWICC=>GWFSWIDAT(Igrid)%SWICC
        SWICUMCR=>GWFSWIDAT(Igrid)%SWICUMCR
        SWICUMCC=>GWFSWIDAT(Igrid)%SWICUMCC
        NUTOP=>GWFSWIDAT(Igrid)%NUTOP
        NUBOT=>GWFSWIDAT(Igrid)%NUBOT
        BRHS=>GWFSWIDAT(Igrid)%BRHS
        QZEXTRA=>GWFSWIDAT(Igrid)%QZEXTRA
        DUM=>GWFSWIDAT(Igrid)%DUM
        RHSFRESH=>GWFSWIDAT(Igrid)%RHSFRESH
        HCOFOLD=>GWFSWIDAT(Igrid)%HCOFOLD
        IPLPOS=>GWFSWIDAT(Igrid)%IPLPOS
        IZONENR=>GWFSWIDAT(Igrid)%IZONENR
        SWIDE4=>GWFSWIDAT(Igrid)%SWIDE4
        SWIPCG=>GWFSWIDAT(Igrid)%SWIPCG
C
      RETURN
      END SUBROUTINE SGWF2SWI1PNT
C
C-------SAVE POINTERS TO SWI DATA FOR A GRID
      SUBROUTINE SGWF2SWI1PSV(Igrid)
      USE GWFSWIMODULE
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: Igrid
C
        GWFSWIDAT(Igrid)%NPLN=>NPLN
        GWFSWIDAT(Igrid)%ISTRAT=>ISTRAT
        GWFSWIDAT(Igrid)%ISWIZT=>ISWIZT
        GWFSWIDAT(Igrid)%NPRS=>NPRS
        GWFSWIDAT(Igrid)%NZONES=>NZONES
        GWFSWIDAT(Igrid)%NLAYSWI=>NLAYSWI
!        GWFSWIDAT(Igrid)%NODES=>NODES
!        GWFSWIDAT(Igrid)%NHALFU=>NHALFU
!        GWFSWIDAT(Igrid)%NHALFL=>NHALFL
        GWFSWIDAT(Igrid)%USEDE45=>USEDE45
!        GWFSWIDAT(Igrid)%ID4DIM=>ID4DIM
!        GWFSWIDAT(Igrid)%ID4DIR=>ID4DIR
!        GWFSWIDAT(Igrid)%NBWGRD=>NBWGRD
!        GWFSWIDAT(Igrid)%ITMX=>ITMX
!        GWFSWIDAT(Igrid)%MXITER=>MXITER
!        GWFSWIDAT(Igrid)%MXUP=>MXUP
!        GWFSWIDAT(Igrid)%MXLOW=>MXLOW
!        GWFSWIDAT(Igrid)%MXBW=>MXBW
!        GWFSWIDAT(Igrid)%MXEQ=>MXEQ
!        GWFSWIDAT(Igrid)%MUTD4=>MUTD4
!        GWFSWIDAT(Igrid)%NBWL=>NBWL
!        GWFSWIDAT(Igrid)%NUPL=>NUPL
!        GWFSWIDAT(Igrid)%NLOWL=>NLOWL
!        GWFSWIDAT(Igrid)%NLOW=>NLOW
!        GWFSWIDAT(Igrid)%NEQ=>NEQ
!        GWFSWIDAT(Igrid)%NUP=>NUP
!        GWFSWIDAT(Igrid)%NBW=>NBW
!        GWFSWIDAT(Igrid)%DELTL=>DELTL
        GWFSWIDAT(Igrid)%TOESLOPE=>TOESLOPE
        GWFSWIDAT(Igrid)%TIPSLOPE=>TIPSLOPE
        GWFSWIDAT(Igrid)%ZETAMIN=>ZETAMIN
        GWFSWIDAT(Igrid)%DELZETA=>DELZETA
!        GWFSWIDAT(Igrid)%AU=>AU
!        GWFSWIDAT(Igrid)%IUPP=>IUPP
!        GWFSWIDAT(Igrid)%AL=>AL
!        GWFSWIDAT(Igrid)%IEQP=>IEQP
!        GWFSWIDAT(Igrid)%SWID4B=>SWID4B
!        GWFSWIDAT(Igrid)%SWILRCH=>SWILRCH
!        GWFSWIDAT(Igrid)%SWIHDCG=>SWIHDCG
        GWFSWIDAT(Igrid)%IBO=>IBO
        GWFSWIDAT(Igrid)%SWIHCOF=>SWIHCOF
        GWFSWIDAT(Igrid)%SWISOLVCR=>SWISOLVCR
        GWFSWIDAT(Igrid)%SWISOLVCC=>SWISOLVCC
        GWFSWIDAT(Igrid)%SWISOLVCV=>SWISOLVCV
        GWFSWIDAT(Igrid)%ZETA=>ZETA
        GWFSWIDAT(Igrid)%ZETAOLD=>ZETAOLD
        GWFSWIDAT(Igrid)%SSZ=>SSZ
        GWFSWIDAT(Igrid)%EPS=>EPS
        GWFSWIDAT(Igrid)%NUS=>NUS
        GWFSWIDAT(Igrid)%DELNUS=>DELNUS
        GWFSWIDAT(Igrid)%NUPLANE=>NUPLANE
        GWFSWIDAT(Igrid)%SWICR=>SWICR
        GWFSWIDAT(Igrid)%SWICC=>SWICC
        GWFSWIDAT(Igrid)%SWICUMCR=>SWICUMCR
        GWFSWIDAT(Igrid)%SWICUMCC=>SWICUMCC
        GWFSWIDAT(Igrid)%NUTOP=>NUTOP
        GWFSWIDAT(Igrid)%NUBOT=>NUBOT
        GWFSWIDAT(Igrid)%BRHS=>BRHS
        GWFSWIDAT(Igrid)%QZEXTRA=>QZEXTRA
        GWFSWIDAT(Igrid)%DUM=>DUM
        GWFSWIDAT(Igrid)%RHSFRESH=>RHSFRESH
        GWFSWIDAT(Igrid)%HCOFOLD=>HCOFOLD
        GWFSWIDAT(Igrid)%IPLPOS=>IPLPOS
        GWFSWIDAT(Igrid)%IZONENR=>IZONENR
        GWFSWIDAT(Igrid)%SWIDE4=>SWIDE4
        GWFSWIDAT(Igrid)%SWIPCG=>SWIPCG
C
      RETURN
      END SUBROUTINE SGWF2SWI1PSV
